* hiroshima_reducedform.do

graph set window fontface "Arial"

use "../../data/processed data/hiroshima_main_1933_75.dta", clear

global tabledir "../../output/table" 

global figuredir "../../output/figure"


ssc install psacalc, replace 


* *******************************************
* Summary stats (Table A.2, A.3)
* *******************************************

* Small-damaged + no damage buildings
gen littledamage_build = smalldamage_build + nodamage_build

* Additional controls 
*(land use regulation, old city dummy, distance to river, administrative boundary (city vs prefecture))
* Dummy of whether the block includes the block (within 10 meters)
rename distance dist_ri
gen river_dummy = (dist_ri<10)
* Dummies for land-use regulation. We assume that manufacturing, unspecified, outsidezone are not regulated (see メモ in 戦前用途地域規制フォルダ) -> Control for housing, commercial, rosen commercial. Baseline=free areas
global z_additional housing commercial rosencomme oldcity kukaku river_dummy


* full sample
estpost summarize pop36 popdens36 pop45bomb popdens45bomb pop49 popdens49 pop50 popdens50 pop51 popdens51 pop60 popdens60 pop75 popdens75 emp38 empdens38 emp45bomb empdens45bomb emp53 empdens53 emp66 empdens66 emp75 empdens75  lnQ_50 lnQ_55 lnQ_60 lnQ_65 lnQ_70 lnQ_75 size cbd_dist 
esttab using "$tabledir/summary_stat_basicdata.tex", label replace cells("mean(fmt(%8.3g)) sd(fmt(%8.3g)) count(fmt(%5.0g))") nomtitle nonumber

* Summary statistics for fundamentals. Table A.3
estpost summarize size cbd_dist port_dist dis_sta dis_wat dis_cass altitude slope river_dummy badsoil lat lon pretrend tdest_build halfdamage_build smalldamage_build nodamage_build housing commercial rosencomme kukaku oldcity 
esttab using "$tabledir/summary_stat_fundamentals.tex", label replace cells("mean(fmt(%8.3g)) sd(fmt(%8.3g)) count(fmt(%5.0g))") nomtitle nonumber

* Within 3km
estpost summarize pop36 popdens36 pop45bomb popdens45bomb pop49 popdens49 pop50 popdens50 pop51 popdens51 pop60 popdens60 pop75 popdens75 emp38 empdens38 emp45bomb empdens45bomb emp53 empdens53 emp66 empdens66 emp75 empdens75  lnQ_50 lnQ_55 lnQ_60 lnQ_65 lnQ_70 lnQ_75 size cbd_dist   if cbd_dist<3.0 
esttab using "$tabledir/summary_stat_basicdata_3km.tex", label replace cells("mean(fmt(%8.3g)) sd(fmt(%8.3g)) count(fmt(%5.0g))") nomtitle nonumber

* Summary statistics for fundamentals. Table A.3
estpost summarize size cbd_dist port_dist dis_sta dis_wat dis_cass altitude slope river_dummy badsoil lat lon pretrend tdest_build halfdamage_build smalldamage_build nodamage_build housing commercial rosencomme kukaku oldcity  if cbd_dist<3.0 
esttab using "$tabledir/summary_stat_fundamentals_3km.tex", label replace cells("mean(fmt(%8.3g)) sd(fmt(%8.3g)) count(fmt(%5.0g))") nomtitle nonumber

* ****************************************************************************
* ************ Figures relating destruction level to CBD dist  ***************
* ****************************************************************************
* Figure A.2(a)
#d;
twoway (scatter tdest_build cbd_dist, sort color(black) msize(small)),
	legend(off)
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.7)
	xtitle("Distance to CBD (km)", size(6)) xlabel(, angle(horizontal) labsize(large))
	ytitle("Percentage of destroyed buildings by bombing", size(6)) ylabel(, angle(horizontal) labsize(large));
	graph export "$figuredir/Destroybuild_CBDdistance.pdf", as(pdf)  replace ;
#d cr

* Figure A.2(b)
#d;
twoway (scatter death_rate cbd_dist, sort color(black) msize(small)),
	legend(off)
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.7)
	xtitle("Distance to CBD (km)", size(6)) xlabel(, angle(horizontal) labsize(large))
	ytitle("Percentage of death by bombing", size(6)) ylabel(, angle(horizontal) labsize(large));
	graph export "$figuredir/Destroypeople_CBDdistance.pdf", as(pdf)  replace ;
#d cr

* ***************************************************************************
* Figure -- population density in 1936-60 (Adjusting the total population) **
* ***************************************************************************
egen tot_pop36 = total(pop36)
gen pop36_adj = pop36*100000/tot_pop36
egen tot_pop45 = total(pop45bomb)
gen pop45_adj = pop45bomb*100000/tot_pop45
egen tot_pop50 = total(pop50)
gen pop50_adj = pop50*100000/tot_pop50
egen tot_pop60 = total(pop60)
gen pop60_adj = pop60*100000/tot_pop60

foreach i in 36 45 50 60{
	gen popdens`i'_adj = pop`i'_adj/size
	gen l_popdens`i'_adj = ln(pop`i'_adj/size)
}


**  1950 population, accounting for pre-trend in the pre-war period.
egen tot_pop50_prewartrend = total(pop50_prewartrend)
gen pop50_prewartrend_adj = pop50_prewartrend*100000/tot_pop50_prewartrend
gen popdens50_prewartrend_adj = pop50_prewartrend_adj/size
gen l_popdens50_prewartrend_adj = log(popdens50_prewartrend_adj)

* Share version
gen pop36_sh = pop36/tot_pop36
gen pop45_sh = pop45bomb/tot_pop45
gen pop50_sh = pop50/tot_pop50
gen pop60_sh = pop60/tot_pop60
foreach i in 36 45 50 60{
	gen popdens`i'_sh = pop`i'_sh/size
	gen l_popdens`i'_sh = log(pop`i'_sh/size)
}
gen pop50_prewartrend_sh = pop50_prewartrend/tot_pop50_prewartrend
gen popdens50_prewartrend_sh = pop50_prewartrend_sh/size
gen l_popdens50_prewartrend_sh = log(popdens50_prewartrend_sh)


** Figure 3 in paper
#d;
twoway (lpoly l_popdens36_adj cbd_dist [aw=size], bwidth(0.4) sort color(black) lwidth(thick) lpattern(shortdash))
(lpoly l_popdens45_adj cbd_dist [aw=size], bwidth(0.4) sort color(cranberry) lwidth(thick) lpattern(dash)) 
(lpoly l_popdens50_adj cbd_dist [aw=size], bwidth(0.4) sort color(black) lwidth(thick) lpattern(line)) 
(lpoly l_popdens50_prewartrend_adj cbd_dist [aw=size], bwidth(0.4) sort color(black) lwidth(thick) lpattern(dash_dot)),
	legend(order(1 "1936" 2 "1945 (after bombing)" 3 "1950" 4 "Predicted 1950 (following pre-trend without bombing)") row(2) size(large))
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.7)
	xtitle("Distance to CBD (km)", size(6.0)) xlabel(, angle(horizontal) labsize(large))
	ytitle("Log population density (per square km)", size(6.0)) ylabel(, angle(horizontal) labsize(large));
	graph export "$figuredir/LogPopdensityPREWARTREND_CBDdistance_diffyears_lpoly3.pdf", as(pdf)  replace ;
#d cr


* *****************************************
* Figures -- employment density 
* *****************************************


egen tot_emp38 = total(emp38)
gen emp38_adj = emp38*100000/tot_emp38
egen tot_emp45 = total(emp45bomb)
gen emp45_adj = emp45bomb*100000/tot_emp45
egen tot_emp53 = total(emp53)
gen emp53_adj = emp53*100000/tot_emp53
egen tot_emp66 = total(emp66)
gen emp66_adj = emp66*100000/tot_emp66

foreach i in 38 45 53 66{
	gen empdens`i'_adj = emp`i'_adj/size
	gen l_empdens`i'_adj = log(emp`i'_adj/size)
}


* Log employment density (same bandwidth as population graphs)
* Figure A.7 in paper
#d;
twoway (lpoly l_empdens38_adj cbd_dist [aw=size], bwidth(0.34) sort color(black) lwidth(thick) lpattern(shortdash)) 
(lpoly l_empdens45_adj cbd_dist [aw=size], bwidth(0.34) sort color(cranberry) lwidth(thick) lpattern(dash))
(lpoly l_empdens53_adj cbd_dist [aw=size], bwidth(0.34) sort color(black) lwidth(thick) lpattern(dash_dot)),
	legend(order(1 "1938" 2 "1945 (after bombing)" 3 "1950") row(1) size(large))
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.7)
	xtitle("Distance to CBD (km)", size(6.0))  xlabel(, angle(horizontal) labsize(large))
	ytitle("Log Employment density (per square km)", size(6.0)) ylabel(, angle(horizontal) labsize(large));
	graph export "$figuredir/LogEmploymentADJUSTED_CBDdistance_diffyears_lpolysmoother.pdf", as(pdf)  replace ;
#d cr


gen dpopdens65 = l_popdens65 - l_popdens45bomb
gen dpopdens60 = l_popdens60 - l_popdens45bomb
gen dpopdens51 = l_popdens51 - l_popdens45bomb
gen dpopdens50 = l_popdens50 - l_popdens45bomb
gen dpopdens49 = l_popdens49 - l_popdens45bomb
gen dpopdens70 = l_popdens70 - l_popdens45bomb
gen dpopdens75 = l_popdens75 - l_popdens45bomb
gen dpopdens45 = l_popdens45bomb- l_popdens36

gen dpopdens5055 = l_popdens55 - l_popdens50
gen dpopdens6055 = l_popdens60 - l_popdens55
gen dpopdens6560 = l_popdens65 - l_popdens60
gen dpopdens7065 = l_popdens70- l_popdens65




* ***************************************************************************************************************************
* ***************************************************************************************************************************
* **********************************  Main specification for Population   ***************************************************
* ***************************************************************************************************************************
* ***************************************************************************************************************************

** To align with other specification, 
replace cbd_dist = 1000*cbd_dist
replace port_dist = 1000*port_dist

gen ln_dsta = ln(dis_sta)
gen ln_dwat = ln(dis_wat)
gen ln_dcass = ln(dis_cass)
gen ln_dport = ln(port_dist)
gen ln_altitude = ln(altitude)

gen ln_cbddist = ln(cbd_dist)
gen cbd_dist2 = cbd_dist^2
gen ln_epidist = ln(epi_dist)

* **********************************
* ************ Controls ************
* **********************************

* Log distance versions of control
* 1st-nature controls (natural fundamentals): 
global z_1st dis_wat ln_altitude slope slope2 badsoil lat lon latlon river_dummy
global z_1st_ln ln_dwat ln_altitude slope slope2 badsoil lat lon latlon river_dummy

* 2nd-nature controls (manmade fundamentals):
global z_2nd dis_sta dis_cass  port_dist  littledamage_build housing commercial rosencomme oldcity kukaku
global z_2nd_ln ln_dsta ln_dcass ln_dport littledamage_build housing commercial rosencomme oldcity kukaku

* Total  
global z dis_wat ln_altitude slope slope2 badsoil lat lon latlon dis_sta dis_cass  port_dist  littledamage_build housing commercial rosencomme oldcity kukaku river_dummy
global z_ln ln_dwat ln_altitude slope slope2 badsoil lat lon latlon ln_dsta ln_dcass ln_dport littledamage_build housing commercial rosencomme oldcity kukaku river_dummy
global z_trend pretrend pretrend_q

global z_nolatlon ln_dwat ln_altitude slope slope2 badsoil ln_dsta ln_dcass ln_dport littledamage_build housing commercial rosencomme oldcity kukaku river_dummy
global z_1st_nolatlon ln_dwat ln_altitude slope slope2 badsoil river_dummy

global z_nolatlon_nobuild ln_dwat ln_altitude slope slope2 badsoil ln_dsta ln_dcass ln_dport river_dummy

********************************************************************************************
* Figures for Population Density Changes over Different Time Horizons
* The fitted line is an unweighted simple OLS regression.
********************************************************************************************
* Figure B.1(a)
#d;
twoway (scatter dpopdens50 dpopdens45  [w=popdens36], msymbol(circle) color(black%15)) 
(lfit dpopdens50 dpopdens45, sort color(cranberry) lwidth(thick))
(function y=-x, range(dpopdens45) color(black) lwidth(thick) lpattern(dash)),
    legend(order(2 "Fitted line" 3 "Slope of -1") row(1) size(large))
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.7)
	xtitle("Change in log population density, 1936-45 (after bombing)", size(6)) xlabel(, labsize(large))
	ytitle("Change in log population density, 1945-50", size(6)) ylabel(-1 0 1 2 3 4 5, angle(horizontal) labsize(large));
	graph export "$figuredir/Popdensity_Change_DW_1stperiod.pdf", as(pdf)   replace ;
#d cr

* Figure B.1(b)
#d;
twoway (scatter dpopdens60 dpopdens45  [w=popdens36], msymbol(circle) color(black%15)) 
(lfit dpopdens60 dpopdens45, sort color(cranberry) lwidth(thick))
(function y=-x, range(dpopdens45) color(black) lwidth(thick) lpattern(dash)),
    legend(order(2 "Fitted line" 3 "Slope of -1") row(1) size(large))
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.7)
	xtitle("Change in log population density, 1936-45 (after bombing)", size(6)) xlabel(, labsize(large))
	ytitle("Change in log population density, 1945-60", size(6)) ylabel(-1 0 1 2 3 4 5, angle(horizontal) labsize(large));
	graph export "$figuredir/Popdensity_Change_DW_2ndperiod.pdf" , as(pdf)  replace ;
#d cr

* ****************************
** 1950 regressions (Main) ***
* ****************************

reg dpopdens50 dpopdens45, vce(hc3)    /* No Controls */  
est sto dpopdens50_1
test dpopdens45 = -1

predict er1, resid
reg er1 ln_cbddist     /* Error term of the simple regression. Not correlated with log cbd distance */ 

reg dpopdens50 dpopdens45 ${z_1st_ln} , vce(hc3)  /*1st Nature Controls*/ 
est sto dpopdens50_2
test dpopdens45 = -1

reg dpopdens50 dpopdens45 ${z_2nd_ln}, vce(hc3)  /*2nd Nature Controls*/ 
est sto dpopdens50_3
test dpopdens45 = -1

reg dpopdens50 dpopdens45 ${z_1st_ln} ${z_2nd_ln}, vce(hc3)  /* 1st & 2nd Nature Controls */
est sto dpopdens50_4
test dpopdens45 = -1

* Oster's bound. We assume the maximum R^2 is the R^2 under controlled regression plus
* the incremental R^2 from the uncontrolled to controlled one
reg dpopdens50 dpopdens45 ${z_1st_ln} ${z_2nd_ln}, vce(hc3)  /* 1st & 2nd Nature Controls */
psacalc beta dpopdens45, rmax(0.911)

reg dpopdens50 dpopdens45 ${z_1st_ln} ${z_2nd_ln} ${z_trend}, vce(hc3)  /* Add Pre-trend */ 
est sto dpopdens50_5
test dpopdens45 = -1

* Robustness: add cubic trends.
reg dpopdens50 dpopdens45 ${z_1st_ln} ${z_2nd_ln} ${z_trend} pretrend_c, vce(hc3)  /* Add Pre-trend */ 

reg dpopdens50 dpopdens45 ${z_1st_ln} ${z_2nd_ln} if cbd_dist<3000, vce(hc3) /* Focus on 3 km from CBD */ 
est sto dpopdens50_6
test dpopdens45 = -1

* Oster's bound. We assume the maximum R^2 is the R^2 under controlled regression plus
* the incremental R^2 from the uncontrolled to controlled one
reg dpopdens50 dpopdens45 ${z_1st_ln} ${z_2nd_ln} if cbd_dist<3000, vce(hc3) /* Focus on 3 km from CBD */ 
psacalc beta dpopdens45, rmax(0.934)

** Table 1 in paper
esttab dpopdens50_1 dpopdens50_2 dpopdens50_3 dpopdens50_4 dpopdens50_5 dpopdens50_6 ///
using "$tabledir/ols_regtable_popdens.tex", se r2 star(c 0.1 b 0.05 a 0.01) b(4) keep(dpopdens45) replace


* **********************************************
***  OLS with 1936 and 1945 in RHS ***
* Table B.1 in paper
* **********************************************

reg l_popdens50 l_popdens36 l_popdens45bomb, vce(hc3)
est sto dp50_1
lincom l_popdens36 + l_popdens45bomb 

reg l_popdens50 l_popdens36 l_popdens45bomb ${z_1st_ln}, vce(hc3)
est sto dp50_2
lincom l_popdens36 + l_popdens45bomb 

reg l_popdens50 l_popdens36 l_popdens45bomb ${z_2nd_ln}, vce(hc3)
lincom l_popdens36 + l_popdens45bomb 

reg l_popdens50 l_popdens36 l_popdens45bomb ${z_1st_ln} ${z_2nd_ln}, vce(hc3)
est sto dp50_3
lincom l_popdens36 + l_popdens45bomb 


reg l_popdens50 l_popdens36 l_popdens45bomb ${z_1st_ln} ${z_2nd_ln} ${z_trend}, vce(hc3)
lincom l_popdens36 + l_popdens45bomb 

reg l_popdens50 l_popdens36 l_popdens45bomb ${z_1st_ln} ${z_2nd_ln} if cbd_dist<3000, vce(hc3)
est sto dp50_4
lincom l_popdens36 + l_popdens45bomb

* Only fundamentals specifications in the RHS
reg l_popdens50  ${z_1st_ln} ${z_2nd_ln}, vce(hc3)
est sto dp50_5

reg l_popdens50  ${z_1st_ln} ${z_2nd_ln} if cbd_dist<3000, vce(hc3)
est sto dp50_6

esttab dp50_1 dp50_2 dp50_3 dp50_4 dp50_5 dp50_6 using "$tabledir/ols_regtable_popdens_me.tex", se r2 star(c 0.1 b 0.05 a 0.01) b(4) keep(l_popdens36 l_popdens45bomb  ${z_1st_ln} ${z_2nd_ln}) replace


******************************************************************************** 
**********************************************************************************
************************** ROBUSTNESS CHECK  ************************************
*********************************************************************************
*********************************************************************************

*********************************************
*** (1) The effect of public housing   ******
* Table B.4 in paper
*********************************************
replace households_house=0 if households_house==.

* Conditional on having some pub housing, is pub housing located near cbd?
gen l_households_house = log(households_house)
reg households_house cbd_dist, vce(hc3)
reg l_households_house cbd_dist, vce(hc3)
est sto pubhousing_dist1
replace households_house=0 if households_house==.

* Unconditional on having pub house, is it more likely to locate near cbd?
reg households_house cbd_dist, vce(hc3)
poisson households_house cbd_dist, vce(robust)
est sto pubhousing_dist2

reg dpopdens50 dpopdens45  households_house, vce(hc3)
est sto pubhousing_dw1
test dpopdens45 = -1
reg dpopdens50 dpopdens45 ${z_ln} households_house, vce(hc3)
est sto pubhousing_dw2
test dpopdens45 = -1
reg dpopdens50 dpopdens45 ${z_ln} households_house if cbd_dist<3000, vce(hc3)
est sto pubhousing_dw3
test dpopdens45 = -1

* Table B4
esttab pubhousing_dist1 pubhousing_dw1 pubhousing_dw2 pubhousing_dw3 using "$tabledir/pubhousing_dw.tex", se r2 star(c 0.1 b 0.05 a 0.01) b(4) keep(cbd_dist dpopdens45 households_house) replace

********************************************************************************************
** (2) Incorporating the effect of nearby blocks (spatial econometric specification) *******
* Table B.3 in paper
********************************************************************************************
ssc install spgen, replace
* Create the spatial lag of the main explanatory variable and controls. Spatial decay parameter is set to 4.32. This lag excludes the own block (see spgen documentation).
* Ahlfeldt et al 2015 has delta=0.36 when distance is measured in minutes. Assuming 60min=5km by walking distance, we get 60*0.36=4.32
spgen dpopdens45 ln_dwat ln_dcass ln_dport ln_dsta ln_altitude slope badsoil pretrend littledamage_build lat lon latlon housing commercial rosencomme oldcity kukaku river_dummy, lat(lat) lon(lon) swm(exp 4.32) dist(100) dunit(km)
reg dpopdens50 dpopdens45 splag1_dpopdens45_e, vce(hc3)
test dpopdens45 = -1
reg dpopdens50 dpopdens45 splag1_dpopdens45_e ${z_ln}, vce(hc3)
est sto speconom_dw1
test dpopdens45 = -1

* Control variables defined by the spatial lags.
global z_splag splag1_ln_dport_e splag1_ln_dwat_e splag1_ln_dsta_e splag1_ln_dcass_e splag1_ln_altitude_e splag1_slope_e splag1_slope_e2 splag1_badsoil_e splag1_littledamage_build_e splag1_lat_e splag1_lon_e splag1_latlon_e splag1_housing_e splag1_commercial_e splag1_rosencomme_e splag1_oldcity_e splag1_kukaku_e splag1_river_dummy_e
gen splag1_slope_e2 = splag1_slope_e^2
gen splag1_pretrend_e2 = splag1_pretrend_e^2


* With controls, the baseline characteristics and spatial lag
reg dpopdens50 dpopdens45 ${z_splag}, vce(hc3)
test dpopdens45 = -1
reg dpopdens50 dpopdens45 ${z_ln} ${z_splag}, vce(hc3)
est sto speconom_dw2
test dpopdens45 = -1
reg dpopdens50 dpopdens45 splag1_dpopdens45_e ${z_ln} ${z_splag}, vce(hc3)

* Spatial lag of 1945 population (market access term?)
* Addresses the possibility that the center is accessible to surrounding areas.
* Delta=0.05 is based on the delta term in the previous version (January 2023 version. rho/sigma=8, delta*rho/sigma=0.4)
spgen pop45bomb, lat(lat) lon(lon) swm(exp 0.05) dist(100) dunit(km) nostd
reg dpopdens50 dpopdens45 splag1_pop45bomb, vce(hc3)
test dpopdens45 = -1
reg dpopdens50 dpopdens45 splag1_pop45bomb ${z_ln}, vce(hc3)
est sto speconom_dw3
test dpopdens45 = -1
* Considering own population size hardly afffects the result.
gen splag_alt_pop45bomb = splag1_pop45bomb + pop45bomb
reg dpopdens50 dpopdens45 splag_alt_pop45bomb ${z_ln}, vce(hc3)
* Normalizing the density
drop splag1_pop45bomb
spgen pop45bomb, lat(lat) lon(lon) swm(exp 0.05) dist(100) dunit(km)
reg dpopdens50 dpopdens45 splag1_pop45bomb ${z_ln}, vce(hc3)

reg dpopdens50 dpopdens45 splag1_dpopdens45_e splag1_pop45bomb ${z_ln} ${z_splag}, vce(hc3)


* Adding spatial lag of area size (a measure of geographic centrality)
spgen size, lat(lat) lon(lon) swm(exp 0.05) dist(100) dunit(km) nostd
reg dpopdens50 dpopdens45 splag1_size_e , vce(hc3)
reg dpopdens50 dpopdens45 splag1_size_e splag1_pop45bomb_e, vce(hc3)
reg dpopdens50 dpopdens45 splag1_size_e ${z_ln}, vce(hc3)
est sto speconom_dw4
test dpopdens45 = -1
reg dpopdens50 dpopdens45 splag1_size_e splag1_pop45bomb_e  ${z_ln}, vce(hc3)
reg dpopdens50 dpopdens45 splag1_size_e splag1_pop45bomb_e  ${z_splag}, vce(hc3)

reg dpopdens50 dpopdens45 splag1_dpopdens45_e splag1_pop45bomb splag1_size_e ${z_splag} ${z_ln}, vce(hc3)
est sto speconom_dw5
test dpopdens45 = -1


esttab speconom_dw1 speconom_dw2 speconom_dw3 speconom_dw4 speconom_dw5 using "$tabledir/speconom_dw.tex", se r2 star(c 0.1 b 0.05 a 0.01) b(4) keep(dpopdens45) replace


* *************************
* *** (3) Employment ****** 
* Table B.2 in paper (1953 employment is assumed to be proportional to 1950 employment)
* *************************
gen dempdens66 = l_empdens66 - l_empdens45bomb 
gen dempdens57 = l_empdens57 - l_empdens45bomb 
gen dempdens53 = l_empdens53 - l_empdens45bomb 
gen dempdens45 = l_empdens45bomb - l_empdens38

* ols (for 1953 employment. Same result for 1950 employment)
reg dempdens53 dempdens45, vce(hc3)
est sto dempdens53_1
test dempdens45 = -1
reg dempdens53 dempdens45 ${z_1st_ln}, vce(hc3)
est sto dempdens53_2
test dempdens45 = -1
reg dempdens53 dempdens45 ${z_2nd_ln}, vce(hc3)
est sto dempdens53_3
test dempdens45 = -1
reg dempdens53 dempdens45 ${z_1st_ln} ${z_2nd_ln}, vce(hc3)
est sto dempdens53_4
test dempdens45 = -1
reg dempdens53 dempdens45 ${z_1st_ln} ${z_2nd_ln} ${z_trend}, vce(hc3)
est sto dempdens53_5
test dempdens45 = -1
reg dempdens53 dempdens45 ${z_1st_ln} ${z_2nd_ln} if cbd_dist<3000, vce(hc3)
est sto dempdens53_6
test dempdens45 = -1

esttab dempdens53_1 dempdens53_2 dempdens53_3 dempdens53_4 dempdens53_5 dempdens53_6 using "$tabledir/ols_regtable_empdens53.tex", se r2 star(c 0.1 b 0.05 a 0.01) b(4) keep(dempdens45) replace

* Correlation between establishment density and employment density in 1966
* Figure A.8(a)
#d;
twoway (scatter l_empdens66 l_estdens66, sort color(black) msize(small)) 
(lfit l_empdens66 l_estdens66, sort color(cranberry) lwidth(medthick)),
	scheme(s1color) plotregion(lwidth(none) ilwidth(none))
    legend(off)
	xtitle("Log establishment density in 1966")
	ytitle("Log employment density in 1966") ylabel(, angle(horizontal));
	graph export "$figuredir/est_emp_1966.pdf", as(pdf)   replace ;
#d cr


* ******************************************************************************
* ************ FIGURE FOR POPULATION AND EMPLOYMENT CHANGE *********************
* ******************************************************************************
* Population 
* Figure 2(a) in paper
rename pop45bomb pop45 
reshape long pop, i(block) j(year)

* To calculate share withim 1km
gen pop2 = pop if cbd_dist < 1000

collapse (sum) pop pop2, by(year)

replace pop = pop/1000
replace pop2 = pop2/1000
gen popcshare = pop2/pop

tsset year 

#d;
twoway (tsline pop, xtick(33 34 35 36 45 49 50 51 52 53 55 60 65 70 75) xlabel(33 36 45 49 51 53 55 60 65 70 75)  ylabel(, format(%8.0g)) recast(connected) color(black) lwidth(thick))
(tsline pop2, xtick(33 34 35 36 45 49 50 51 52 53 55 60 65 70 75) xlabel(33 36 45 49 51 53 55 60 65 70 75)  ylabel(, format(%8.0g)) recast(connected) color(black) lwidth(thick) lpattern(shortdash))
(tsline popcshare, yaxis(2) xtick(33 34 35 36 45 49 50 51 52 53 55 60 65 70 75) xlabel(33 36 45 49 51 53 55 60 65 70 75)  ylabel(, format(%8.0g)) recast(connected) color(cranberry) lwidth(thick) lpattern(dash_dot)),
	legend(order(1 "Total Number" 2 "Number within 1km from CBD" 3 "Share within 1km from CBD (right axis)") row(2) size(3.5))
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.7)
	xtitle("Year", size(6.0)) xlabel(, labsize(medlarge))
	ytitle("Population in Hiroshima City (thousands)", axis(1) size(6.0)) ylabel(, axis(1) angle(horizontal) labsize(large))
	ytitle("Share in 1km from the CBD", color(cranberry) axis(2) size(6.0)) ylabel(, labcolor(cranberry) axis(2) angle(horizontal) labsize(large));
	graph export "$figuredir/Hiroshima_pop_year.pdf" , as(pdf)  replace ;
#d cr

* Employment 
* Figure 2(b) in paper
use "../../data/processed data/hiroshima_main_1933_75.dta", clear
rename emp45bomb emp45 

reshape long emp, i(block) j(year)


gen emp2=emp if cbd_dist<1.0 

collapse(sum) emp emp2, by(year) 
replace emp = emp/1000 
replace emp2 = emp2/1000 
gen empcshare = emp2/emp 

tsset year 

#d;
twoway (tsline emp, xtick(38 45 53 57 63 66 75) xlabel(38 45 53 57 63 66 75)  ylabel(, format(%8.0g)) recast(connected) color(black) lwidth(thick))
(tsline emp2, xtick(38 45 53 57 63 66 75) xlabel(38 45 53 57 63 66 75)  ylabel(, format(%8.0g)) recast(connected) color(black) lwidth(thick) lpattern(shortdash))
(tsline empcshare, yaxis(2) xtick(38 45 53 57 63 66 75) xlabel(38 45 53 57 63 66 75)  ylabel(, format(%8.0g)) recast(connected) color(cranberry) lwidth(thick) lpattern(dash_dot)),
	legend(order(1 "Total Number" 2 "Numbers within 1km from CBD"  3 "Share within 1km from CBD (right axis)") row(2))
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.7)
	xtitle("Year", size(6.0)) xlabel(, labsize(medlarge))
	ytitle("Employment in Hiroshima City (thousands)", axis(1) size(6.0)) ylabel(, axis(1) angle(horizontal) labsize(large))
	ytitle("Share in 1km from CBD", axis(2) size(6.0) color(cranberry)) ylabel(, axis(2) angle(horizontal) labsize(large) labcolor(cranberry));
	graph export "$figuredir/Hiroshima_emp_year.pdf" , as(pdf)  replace ;
#d cr



